Direct Numerical Simulations of Electrophoresis of Charged Colloids 
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We propose a numerical method to simulate electrohydrodynamic phenomena in charged colloidal 
dispersions. This method enables us to compute the time evolutions of colloidal particles, ions, and 
host fluids simultaneously by solving Newton, advection-diffusion, and Navier-Stokes equations 
so that the electrohydrodynamic couplings can be fully taken into account. The electrophoretic 
mobilities of charged spherical particles are calculated in several situations. The comparisons with 
approximation theories show quantitative agreements for dilute dispersions without any empirical 
parameters, however, our simulation predicts notable deviations in the case of dense dispersions. 
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Electrohydrodynamic phenomena are of great impor- 
tance in physical, chemical, and biological science, and 
also in several engineering fields y| ■ In the case of elec- 
trophoresis of charged particles for example, the parti- 
cles start to move on the application of external electric 
fields. The electric double layer, i.e. the cloud of coun- 
terions around charged particles, tends to be deformed 
and its distribution becomes anisotropic because of the 
applied external field and also of the friction between ions 
and fluids. The electrophoretic mobility of a single col- 
loidal particle is then determined by the balance between 
the electrostatic driving force and the hydrodynamic fric- 
tional force acting on the particle. In this situation, the 
time evolutions of the colloidal particles, the ions, and 
the host fluids are described by coupled equations of hy- 
drodynamics (Navier-Stokes) and electrostatics (Poisson) 
with proper boundary conditions imposed on the surfaces 
of the colloidal particles. However, the usual numerical 
techniques of partial differential equations are hopeless 
to deal with dynamical evolutions of many-particle sys- 
tems since the moving particle-fluid boundary condition 
must be treated at every discrete time step. 

In late years, efforts to resolve hydrodynamic interac- 
tions in colloidal dispersions attract much attention. Var- 
ious advanced methods have been proposed such as the 
Stokesian Dynamics (SD) Q, the finite element method 
(FEM) 0, the Lattice Boltzmann method (LBM) 0, 
the Stochastic Rotation Dynamics 0, the Fluid Par- 
ticle Dynamics (FPD) ;6], and yet another method 
which treats solid-fluid interaction efficiently @|. Pi- 
oneering approaches have been proposed also to simu- 
late charged colloidal dispersions without hydrodynam- 
ics 0, 0, 03 El 113 • Although extensions have been done 
to take into account t he hydro dynamics by using SD [13( , 
FEM 0, LBM EmElElli, and FPD E5> quanti- 
tatively reliable simulations have not yet been performed 
successfully for many particle dispersions due to the com- 
plexity of the system. 

Recently, we developed a reliable and efficient numeri- 
cal method, called smoothed profile (SP) method 0,H2]j 



to resolve the hydrodynamic interactions acting on solid 
particles immersed in Newtonian host fluids. In the SP 
method, the original sharp boundaries between colloids 
and host fluids are replaced with diffuse interfaces with 
finite thickness £. This simple modification greatly im- 
prove the performance of numerical computations since 
it enables us to use the fixed Cartesian grid even for the 
problems with moving boundary conditions. 

The SP method is not only applicable to the disper- 
sions in Newtonian fluids, but particularly suitable for 
the particle dispersions in complex fluids. It has already 
been applied successfully to liquid crystal colloidal dis- 
persions [2j3, and charged colloidal dispersions [25| . 
Field-particle hybrid simulations were performed, where 
the average direction of the liquid crystal molecules and 
the density of ions were treated as coarse-grained contin- 
uum objects while colloids were treated explicitly as par- 
ticles. The interaction between fields and particles were 
taken through the diffuse interface. The above meth- 
ods for the dispersions in complex fluids arc, however, 
not yet appropriate for simulating dynamical phenom- 
ena since hydrodynamic effects are completely neglected. 
The purpose of the present study is to establish an effi- 
cient and reliable simulation method applicable for elec- 
trohydrodymanic phenomena such as electrophoresis by 
combining our SP methods for hydrodynamic |2lll22j and 
electrostatic [2^ interactions. 

In the present paper, we briefly outline our numeri- 
cal modeling for charged colloidal dispersions and then 
demonstrate the reliability of the combined SP method 
by comparing our numerical results with classical approx- 
imation theories [^E^ll^l^ . Finally, comparisons are 
made for the electrophoretic mobilities of dense disper- 
sions, where the simulation results show notable devia- 
tions from a mean-field type theory according to the cell 
model 

Let us consider N spherical particles with radius a, the 
mass M p , and the inertia tensor I p in a host fluid con- 
sisting of multi-component ions of species a with charges 
Z a e, where e is the unit charge. The local number density 
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of a ion is C a (f, t) at a time t. The total charge on a col- 
loidal particle is Ze and distributed uniformly on its sur- 
face. The velocity field of the host fluid is v(f, t). The po- 
sition, the translational velocity, and the angular velocity 
of the zth particle are Ri, Vi, and fij, respectively. We 
define the overall profile function <j>(r, t) = J^ i=1 4>i(r, t), 
where 4>i £ [0, 1] is the ith particle's profile field which is 
unity in the particle domain \r — Ri\ < a — £/2, zero in 
the fluid domain \r — Ri\ > a + £/2, and have a contin- 
uous diffuse interface within the thin interface domain 
a — £/2 < \r — Ri\ < a + £/2 whose thickness is £ . 
The mathematical definition of <pi is given in Ref. 211. 
We define the spacial distribution of the surface charge 
eq(r) = Ze\Vcj)\/4:ira 2 using (j>, then the local density of 
the total charge is represented smoothly everywhere in 
the system by p e (r) = ^ a Z a eC a + eq. The complete 
dynamics of the system is obtained by solving the follow- 
ing time evolution equations 0, H^] . 

i) The Navier-Stokes equation: 

P (d t +V-V)v= -Vp + ^V-peV^ + ^e^ + cjjfp, (1) 

with incompressible condition V • v — 0, where p is the 
mass density, p is the pressure, rj is the shear viscosity 
of the host fluid, f el = — E ■ f is the external electric 
potential due to the constant electric field E, and <j)f p 
represents the body force arising from the rigidity of the 
particles jl^. The electrostatic potential ^(r) is to be 
determined by solving the Poisson equation eV 2 ^ = — p e 
with the dielectric constant e of the host fluid. 

ii) The Newton's and Euler's equations of motions: 

R\ = V u MpP^Ff + F?, I p .fk = N?, (2) 

where Ff 1 and N? 1 are the hydrodynamic force and 
torque 22] , and Ff is the force arising from the excluded 
volume of particles which prevents colloids from over- 
lapping. Hereafter, soft-core potential of the truncated 
Lennard- Jones potential is adopted for Ff. We include 
the electric driving force due to E in F^ . 
Hi) Advection-diffusion equation: 

dtC* = -V • C^v + f-W ■ ((I - riri) ■ C* a Vp a ), (3) 

where f a is the ionic friction coefficient, J is the unit 
tensor, and n is a unit- vector field defined by ri — 
— V0/|V0|. In our method, the actual density fields of 
ions are defined as C a (f,t) = (1 — 0(r, £))C*(r, t) us- 
ing the auxiliary density field C*(r, t). This definition 
avoids penetration of ions into colloids explicitly without 
using artificial potentials, which requires smaller time in- 
crements. The operator (J — riri) in Eq.® ensures the 
conservation of C a since the no-penetration condition, 
ri ■ \7p a — is directly assigned at the diffuse interface. 
Then the charge neutrality J p e dr = of the total system 
is guaranteed automatically. Based on the density func- 
tional theory [32I l33| , the chemical potential is defined 
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FIG. 1: Relationship between surface charge \Z\e and dimen- 
sionless zeta potential y (a). Our numerical data follows nicely 
on the analytic solution of the nonlinear PB equation |3^l but 
deviates notably from the Debye-Hiickel linearized theory. 
Dimensionless mobility E m of a single particle is plotted in 
(b) as a function of dimensionless zeta potential y. For com- 
parison, results of Smolchowski, Henry, and O'Brien- White 
for k<2 = 0.5 are plotted. The color contours in (c) and (d) 
represent the total ionic charge density ~}2 a eZ a C a around a 
single particle for (c) Z = -100 and (d) Z = -500. The 
electric field is applied in the horizontal (+x) direction. 

as 

fi a = k B TlnC*+Z a e(* + y ex ), (4) 

where kg is the Boltzmann constant and T is the tem- 
perature. If we set v = in Eq. the equilibrium 
(t — » 00) ionic density is given by the Boltzmann Eq. 

C* a = C a exp[-Z a e(y + V ex )/k B T}. (5) 

The bulk salt concentration C a is related to the De- 
bye screening length kT 1 = 1/^JattXb J2 a ^a^a, where 
Xb = e 2 / AireksT is the Bjerrum length which is approx- 
imately 0.72nm in water at 25 °C. 

Simulations have been performed in a three- 
dimensional cubic box with the periodic boundary con- 
ditions. The linear dimension is L/A = 64, where A 
is the lattice spacing chosen as a unit length, which is 
taken related to the Bjerrum length as A — AttXb- Wc 
use the particle radius a — 5 and the thickness parameter 
£ = 2 throughout the present simulations. The host fluid 
contains 1:1 electrolyte composed of monovalent counte- 
rions (a — +) and coions (a = — ). The unit of the energy 
and the electrostatic potential is ksT and ksT/e, respec- 
tively. The later corresponds to 25.7mV at 25 °C. The 
non-dimensional parameter m a = 2eksT f a /3r]e 2 is set to 
be m+ = m_ = 0.184, which corresponds to the value of 
KC1 solution at 25 °C. Our unit of time r = A 2 f + /k B T 
corresponds to 0.44/isec. 

We first consider a single charged particle moving with 
the drift velocity V = (—V, 0, 0) in a constant electric 
field E = [E, 0, 0). The electrophoretic mobility V/E is 
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related to the zeta potential £, which is defined as the 
electrostatic potential on a slipping plane, as 

V/E = feC/v (6) 

when £ is small The Smolchowski equation / = 1 is 
valid in the limiting case na — > oo |27| . while the Hiickel 
equation / = 2/3 is valid in the opposite case na — > |2sj . 
Henry derived an expression / = /#(/ca) which is valid 
for a general value of na |29| . These equations indicate 
that the mobility is proportional to £, however, this rela- 
tion tends to fail for larger £ where the relaxation effect 
due to deformations of electric double layer becomes no- 
table. O'Brien and White proposed an approximation 
theory which is valid also for larger £ [26| . 

We have performed simulations for electrophoresis of 
a single particle in linear response regimes E < 0.15 and 
compared them with the O'Brien- White theory. A con- 
stant uniform electric field E = 0.1, which corresponds 
to 2.85 x 10 3 V/cm, was applied. The terminal V was 
calculated for 50 < -Z < 750 with k" 1 = 10, corre- 
sponding to the salt concentration 11/xmol/l at 25 °C in 
water. We chose v = r\j p = 5, so the Reynolds number 
Re — aV/v remains small. Both in the O'Brien- White 
theory and our simulations, the zeta potential is assumed 
to be the electrostatic potential at the particle surface, 
(, = ^surface- In our simulations, the surface charges were 
chosen as Z = -50, -100, -200, -300, -400, -500, and 
-750, corresponding to y = 0.525, 1.044, 2.035, 2.927, 
3.692, 4.332, and 5.510, respectively. Here the dimcn- 
sionless zeta potential y = e(/ksT is introduced. A rela- 
tionship between the surface charge \Z\e and the dimcn- 
sionless zeta potential y is shown in Fig^a), where our 
numerical results are plotted together with an analytic 
solution of the nonlinear PB equation [34| and the Debye- 
Hiickel linearized solution £ = |Z|e/47ra 2 eK(l + Ka _1 ). 
We see that our numerical results are consistent with 
the nonlinear PB theory. In FigQIb), the dimension- 
less mobility E m = 3er]V/2ekBTE is plotted as a func- 
tion of the dimensionless zeta potential with na = 0.5. 
It is clearly demonstrated that our method reproduced 
the O'Brien- White theory almost perfectly including the 
nonlinear regime y > 2 only within a few percent error. 
We emphasize that such a precise agreement with the 
theory has never been obtained by any simulation meth- 
ods so far proposed. The distributions of charge density 
due to counterions and coions are shown in Fig^c) for 
y = 1.044 and (d) for y = 3.692. One can see that the 
electric double layer is deformed considerably in the non- 
linear regime (d), while it is almost isotropic in the linear 
regime (c). As is mentioned before, the relaxation effect 
due to the deformed double layer causes the nonlinearity 
in E rn . 

Our simulation method is easily applicable to dense 
dispersions consisting of many particles. We thus exam- 
ined the effect of the particle concentration on the elec- 
trophoretic mobility. The linearized theory for a single 
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FIG. 2: Snapshots of the electrophoresis of dense dispersions 
with (a) FCC, (b) BCC, and (c) random particle configura- 
tions. The color map represents the total ionic charge density 
eZ a C a in a plane perpendicular to z axis. The electric 
field is applied in +x direction normal to (1,0,0) face for FCC 
and BCC. See movies |35|l . 



particle Eq. JSJl is still valid for dense dispersions when E 
is small, however, / is now depending both on na and (f. 
Simulations were carried out with Z = —100 and E = 0.1 
for various particle volume fractions ip = 4ira 3 N/3L 3 to 
calculate f(K,a,ip) = r)V/e£E. The Debye length n^ 1 is 
taken to be 5 and 10 which correspond to na — 1 and 
0.5, respectively. The corresponding salt concentration 
is 44^tmol/l for n^ 1 = 5 and 11/xmol/l for n^ 1 = 10, 
respectively. Figure [21 shows typical snapshots of the sys- 
tems with (a) FCC, (b) BCC, and (c) random config- 
urations |35| . The horizontal color map represents the 
charge density for na = 1 on a cross section perpendicu- 
lar to z axis. In the cases of FCC and BCC, E was ap- 
plied perpendicular to the (1,0,0) and (1,1,1) faces, but 
we obtained very small differences only within 1%. 

The mobility coefficient f{na,p) for na = 1 and 0.5 is 
plotted as a function of ip in Fig. |3| (a) and (b), respec- 
tively. We found that / decreases rapidly with increas- 
ing ip. Furthermore, the overall behavior looks almost 
independent of the particle configurations. A theoreti- 
cal model has been proposed by Levine and Neale for 
the electrophoretic mobility of dense dispersions by us- 
ing the cell model [3(| ■ They assumed the situation that 
a spherical particle with radius a is located at the center 
of a spherical container (cell) with radius b and calcu- 
lated V as a function of ko and tp = (a/6) 3 . Ohshima 
proposed a simpler expression for the mobility coefficient 
/ according to the cell model [U- Ohshima's predic- 
tion is shown in Fig. |3 (a) and (b) together with our 
numerical results. The overall agreement between our 
simulation and Ohshima's theory is better in (a) with a 
smaller Debye length = 5 = a than in (b) with a 
larger one k -1 = 10. In both (a) and (b), the simula- 
tion results tend to be larger than the the theory as tp 
increases. We expect that the deviation is attributable to 
the occurrence of overlapping of the electric double layers 
for larger <p because such an effect is totally neglected in 
the theory. To this end, we estimated the effective radius 
a + kT 1 of the ionically dressed particles and defined the 
effective volume fraction p e ff = 47r(a + k~ 1 ) 3 N/3L 3 = 
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FIG. 3: The mobility coefficient f(na, ip) as a function of 
the volume fraction ip in (a) Ka — 1 and (b) na = 0.5. The 
solid lines rep resent the approximation theory proposed by 
Ohshima |3ll| . The theory is confirmed to be accurate for 
<feff < lj however, tends to deviate from our numerical re- 
sults for ifieff > 1 where overlapping of the electric double 
layers becomes notable. 



(1 + (Ka)~ 1 ) 3 >p. As is clearly seen in Fig|3] (a) and (b), 
our results agree well with Ohshima' theory for ip e f / < 1 
where the effect of overlapping is small. However, for 
feff > 1 where the overlapping of the electric double 
layers becomes large, deviations between our simulations 
and the theory become notable. We emphasize that the 
present study is the first successful simulations which 
provide quantitative data necessary to examine the re- 
liability of the Ohshima's cell model calculations includ- 
ing their boundary conditions for electrophoresis in dense 
colloidal dispersions. Our results are consistent with re- 
cent studies which also devoted to take into account the 
effects of double layer overlapping 0, |3(| . 

In summary, we have developed a unique numerical 
method for simulating electrohydrodynamic phenomena 
in colloidal dispersions. The method was first applied 
to simulate electrophoresis of a single spherical particle, 
and we found that our method can reproduce the re- 
liable analytical theory proposed by O'Brien and White 
quantitatively. Simulations were then performed for elec- 
trophoresis of colloids in dense dispersions, and we com- 
pared them with the theoretical analysis based on the 
cell model. We found that the cell model is reliable when 
overlapping of electric double layers is small but less re- 
liable as overlapping becomes larger. 
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